clear

%set a name for the estimation output
estimation_name='abc_rn';

%%%%%%%%%%%%%%%%%%%%%%%%
%estimation settings
%%%%%%%%%%%%%%%%%%%%%%%%

%set risk neutral or not
risk_neutral=1;         %if set to zero, takes the individual lottery choices for crra utility
                        %if set to one, assumes risk neutrality
                        %if set to two, estimates rho as a free parameter
                        
%set constraints
alpha_zero=0;           %sets alpha to zero
beta_zero=0;            %sets beta to zero
kappa_zero=0;           %sets kappa to zero
delta_zero=1;           %sets delta to zero
gamma_zero=1;           %sets gamma to zero
pms_zero = [alpha_zero beta_zero kappa_zero delta_zero gamma_zero];

alpha_negative=0;       %choose -1 to set alpha positive, 1 to set alpha negative, 0 for no constraint on the sign
beta_negative=0;        %choose -1 to set beta positive, 1 to set beta negative, 0 for no constraint on the sign
kappa_negative=0;       %choose -1 to set kappa positive, 1 to set kappa negative, 0 for no constraint on the sign
delta_negative=0;       %choose -1 to set delta positive, 1 to set delta negative, 0 for no constraint on the sign
gamma_negative=0;       %choose -1 to set gamma positive, 1 to set gamma negative, 0 for no constraint on the sign
lambda_negative=-1;     %choose -1 to set lambda positive, 1 to set lambda negative, 0 for no constraint on the sign
pms_pos = [alpha_negative beta_negative kappa_negative delta_negative gamma_negative lambda_negative];

n_steps = 2;           % number of starting points used per parameter

sts = [risk_neutral pms_zero pms_pos n_steps];

%%%%
% import choices from experiment
filename='mnl_data.csv';
dt=csvread(filename);
file2='risk.csv';       % contains the Eckel-Grossman lottery choices
risk=csvread(file2);

no_trials=18;
n = length(dt)/no_trials; % n is the number of subjects in the data

ind_data = zeros(no_trials,12,n);

parfor i = 1:n

    m1 = no_trials*(i-1)+1; m2 = no_trials*i; id = dt(m1,1); 
      
    subdt = [dt(m1:m2,3:6),dt(m1:m2,7:9),dt(m1:m2,10),dt(m1:m2,11),risk(i)*ones(no_trials,1),id*ones(no_trials,1),dt(m1:m2,2)];
    ind_data(:,:,i)=subdt;
    %ids(i)=id;
end

Results = [];
choice_probs = [];
predictions = [];

parfor i = 1:n 
    for t = 1:no_trials
        hold_out = ind_data(t,:,i);
        est_dt = ind_data(:,:,i);
        est_dt(t,:)=[];

        estimates = hm_indiv_est(est_dt,sts);
        Results = [Results;estimates];
        
        pms=zeros(1,6);
        pms(1:5) = estimates(2:6); %estimates of alpha, beta, kappa, delta, gamma
        pms(6) = estimates(8);     %estimate of lambda 
        pms(7) = estimates(7);     %estimate/value of risk parameter 
        
        probs = hm_probs(pms,hold_out,risk_neutral);
        choice_probs = [choice_probs;probs];
        
        choice = hold_out(8);
        preds = zeros(1,31);
        preds(1:12) = hold_out;
        preds(13:19) = estimates(2:8);
        preds(20) = estimates(9);                          % (LL in case of individual estimations)
        preds(21:28) = probs;
        preds(29) = probs(choice);
        
        [M,I]=max(probs);       
        if I == choice
            correct = 1;
        else 
            correct = 0;
        end
        preds(30) = I;
        preds(31) = correct;
        
        predictions = [predictions;preds];
    end  
end

save(fullfile('output','out_of_sample','mat_files',estimation_name));

T = array2table(predictions);
T_names = {'T','R','P','S','y1','y2','y3','choice','game_type','lot','id','game_no','a_est','b_est','c_est','d_est','e_est','r_est','l_est','ll','p1','p2','p3','p4','p5','p6','p7','p8','p_correct','pred_choice','choice_correct'};
T.Properties.VariableNames(1:31) = T_names;
writetable(T,fullfile('output','out_of_sample','predictions',estimation_name));